﻿#include <learn/alg/alg_covering.h>

using namespace xi;
using namespace xi::alg;

#include <Dense>
#include <Sparse>

namespace
{

template <bool Rat>
std::vector<std::shared_ptr<Geom_BSplineCurve<3, Rat>>> fix_curves_order(
    std::vector<std::shared_ptr<Geom_BSplineCurve<3, Rat>>> &&curves)
{
    std::vector<std::shared_ptr<Geom_BSplineCurve<3, Rat>>> result;
    result.push_back(curves.front());
    curves.erase(curves.begin());

    while (!curves.empty())
    {
        auto it = curves.begin();
        auto mini = it;
        auto distance = std::numeric_limits<double>::max();

        auto c = result.back();
        auto cf = c->first_param();
        auto cl = c->last_param();

        while (it != curves.end())
        {
            auto d = *it;
            auto df = d->first_param();
            auto dl = d->last_param();

            auto d1 = c->value(cl).distance(d->value(df));
            auto d2 = c->value(cl).distance(d->value(dl));

            if (d1 < distance)
            {
                mini = it;
                distance = d1;
            }

            if (d2 < distance)
            {
                mini = it;
                distance = d2;
                d->reverse();
            }

            it++;
        }

        (*mini)->set_pole(0, c->value(cl));

        result.push_back(*mini);
        curves.erase(mini);
    }

    return result;
}

template <bool Rat>
std::vector<std::shared_ptr<Geom_BSplineCurve<3, Rat>>> compat_curves(
    std::vector<std::shared_ptr<Geom_BSplineCurve<3, Rat>>> &&curves, bool elevation = false)
{
    using Vec = VecT<3, double>;
    using VecR = VecT<3 + Rat, double>;

    auto cn = curves.size();
    std::vector<std::shared_ptr<Geom_BSplineCurve<3, Rat>>> result(4);
    if (cn == 2)
    {
        result[0] = curves[0];
        result[2] = curves[1];

        // 获得 4 个顶点位置
        auto p1 = result[0]->value(result[0]->first_param());
        auto p2 = result[0]->value(result[0]->last_param());
        auto p3 = result[2]->value(result[2]->first_param());
        auto p4 = result[2]->value(result[2]->last_param());

        result[1] = Geom_BSplineCurve<3, Rat>::linear(p2, p3);
        result[3] = Geom_BSplineCurve<3, Rat>::linear(p4, p1);
    }
    else if (cn == 3)
    {
        result[0] = curves[0];
        result[1] = curves[1];
        result[2] = curves[2];

        // 创建新的边
        auto p1 = result[2]->value(result[2]->last_param());
        auto p2 = result[0]->value(result[0]->first_param());

        result[3] = Geom_BSplineCurve<3, Rat>::linear(p1, p2);
    }
    else
        result = curves;

    // 调整使得对边具有相同节点向量
    result[0]->reverse();
    Alg_Knot<3>::compat(result[0], result[2]);
    result[0]->reverse();

    result[1]->reverse();
    Alg_Knot<3>::compat(result[1], result[3]);
    result[1]->reverse();

    if (elevation)
    {
        for (auto &c : result)
            c->elevate(2);
    }
    return result;
}

} // namespace

std::shared_ptr<Geom_BSplineSurface<3>> Alg_Coons_Covering::cover(
    std::vector<std::shared_ptr<Geom_BSplineCurve<3>>> &&curves)
{
    auto boundaries = compat_curves(fix_curves_order(std::move(curves)));
    if (boundaries.size() == 0 || boundaries.size() > 4)
        return nullptr;

    // 沿着 cur 线性插值构建切向场
    auto buildDeriv = [&](std::shared_ptr<Geom_BSplineCurve<3>> cur, const Vec3 &sv, const Vec3 &ev) {
        auto deriv = Geom_BSplineCurve<3>::linear(sv, ev, cur->first_param(), cur->last_param());
        deriv->elevate(cur->degree());

        std::vector<double> knots(cur->pole_num() - cur->degree() - 1);
        std::copy(cur->knots().begin() + cur->degree() + 1, cur->knots().begin() + cur->pole_num(), knots.begin());
        deriv->refine_knot(std::move(knots));
        return deriv;
    };

    Vec3 V[8];
    for (std::size_t i = 0; i < 4; i++)
    {
        auto [p1, v1] = boundaries[i]->d1(boundaries[i]->first_param());
        auto [p2, v2] = boundaries[i]->d1(boundaries[i]->last_param());
        V[2 * i] = v1;
        V[2 * i + 1] = v2;
    }

    std::shared_ptr<Geom_BSplineCurve<3>> derivs[4];
    derivs[0] = buildDeriv(boundaries[0], -V[7], V[2]);
    derivs[1] = buildDeriv(boundaries[1], V[1], -V[4]);
    derivs[2] = buildDeriv(boundaries[2], V[3], -V[6]);
    derivs[3] = buildDeriv(boundaries[3], -V[5], V[0]);

    // 累计控制点
    std::vector<Vec3> dstPoles(boundaries[0]->pole_num() * boundaries[1]->pole_num());

    // 构造插值点
    auto buildInterior = [&](std::size_t deg, const Vec3 &p1, const Vec3 &p2, const Vec3 &v1, const Vec3 &v2,
                             Vec3 *poles) {
        if (deg > 2)
        {
            poles[0] = p1;
            poles[1] = p1 + v1 / 3;
            poles[2] = p2 - v2 / 3;
            poles[3] = p2;
        }
        else if (deg == 2)
        {
            // 切向延伸后平均插值
            poles[0] = p1;
            poles[1] = ((p1 + v1 / 2) + (p2 - v2 / 2)) / 2;
            poles[2] = p2;
        }
        else
        {
            poles[0] = p1;
            poles[1] = p2;
        }
    };

    // rc 设置按行或列填充
    auto sumPoles = [&](std::shared_ptr<Geom_BSplineCurve<3>> cur, const Vec3 &p1, const Vec3 &p2, const Vec3 &v1,
                        const Vec3 &v2, std::size_t k, char rc) {
        int deg = cur->degree();
        std::vector<Vec3> poles(std::min(deg + 1, 4));
        buildInterior(deg, p1, p2, v1, v2, poles.data());

        std::vector<double> knots(poles.size() * 2);
        std::fill(knots.begin(), knots.begin() + poles.size(), cur->first_param());
        std::fill(knots.begin() + poles.size(), knots.end(), cur->last_param());

        // 根据边界控制点和切向场控制点，计算中间控制点，得到 Bezier 曲线
        auto bez = std::make_shared<Geom_BSplineCurve<3>>(std::move(poles), std::move(knots), std::min(deg, 3));
        bez->elevate(deg);

        knots.resize(cur->pole_num() - cur->degree() - 1);
        std::copy(cur->knots().begin() + cur->degree() + 1, cur->knots().begin() + cur->pole_num(), knots.begin());
        bez->refine_knot(std::move(knots));

        // 填充控制点
        if (rc == 'r')
        {
            for (std::size_t i = 0; i < bez->pole_num(); i++)
                dstPoles[k * boundaries[1]->pole_num() + i] += bez->pole(i);
        }
        else
        {
            for (std::size_t i = 0; i < bez->pole_num(); i++)
                dstPoles[i * boundaries[1]->pole_num() + k] += bez->pole(i);
        }
    };

    // 构建 S1，S2 并累计
    for (std::size_t i = 0; i < boundaries[0]->pole_num(); i++)
        sumPoles(boundaries[1], boundaries[0]->pole(i), boundaries[2]->pole(boundaries[2]->pole_num() - 1 - i),
                 derivs[0]->pole(i), derivs[2]->pole(derivs[2]->pole_num() - 1 - i), i, 'r');
    for (std::size_t i = 0; i < boundaries[1]->pole_num(); i++)
        sumPoles(boundaries[0], boundaries[3]->pole(boundaries[3]->pole_num() - 1 - i), boundaries[1]->pole(i),
                 derivs[3]->pole(derivs[3]->pole_num() - 1 - i), derivs[1]->pole(i), i, 'c');

    // 构建角点插值曲面 T
    auto buildT = [&](std::size_t tnu, std::size_t tnv) {
        Vec3 bpoles[16];
        for (int i = 0; i < 4; i++)
        {
            int j = (i + 3) % 4;
            int k = (i + 1) % 4;

            // 反转标记
            int s = sign(2 - i);
            buildInterior(boundaries[i]->degree(), boundaries[i]->start_point(), boundaries[i]->end_point(),
                          derivs[j]->end_point() * s, derivs[k]->start_point() * s, &bpoles[i * 4]);
        }

        std::vector<Vec3> poles(tnu * tnv);

        // 上下
        for (std::size_t i = 0; i < tnu; i++)
        {
            poles[i * tnv] = bpoles[i];
            poles[(tnu - i) * tnv - 1] = bpoles[8 + i];
        }

        // 左右
        for (std::size_t i = 0; i < tnv; i++)
        {
            poles[(tnu - 1) * tnv + i] = bpoles[4 + i];
            poles[tnv - i - 1] = bpoles[12 + i];
        }

        // 填充内部
        for (std::size_t i = 1; i < tnu - 1; i++)
        {
            for (std::size_t j = 1; j < tnv - 1; j++)
            {
                poles[i * tnv + j] = (poles[i * tnv] * (tnv - j - 1) + poles[(i + 1) * tnv - 1] * j +
                                      poles[j] * (tnu - i - 1) + poles[(tnu - 1) * tnv + j] * i) /
                                     (tnu + tnv - 2);
            }
        }

        std::vector<double> uKnots(tnu * 2);
        std::fill(uKnots.begin(), uKnots.begin() + tnu, boundaries[0]->first_param());
        std::fill(uKnots.begin() + tnu, uKnots.end(), boundaries[0]->last_param());

        std::vector<double> vKnots(tnv * 2);
        std::fill(vKnots.begin(), vKnots.begin() + tnv, boundaries[1]->first_param());
        std::fill(vKnots.begin() + tnv, vKnots.end(), boundaries[1]->last_param());

        return std::make_shared<Geom_BSplineSurface<3>>(std::move(poles), std::move(uKnots), std::move(vKnots), tnu - 1,
                                                        tnv - 1);
    };

    std::size_t tnu = std::min(4, int(boundaries[0]->degree()) + 1);
    std::size_t tnv = std::min(4, int(boundaries[1]->degree()) + 1);
    auto T = buildT(tnu, tnv);

    // 升阶并节点插入
    T->elevate(boundaries[0]->degree(), boundaries[1]->degree());

    std::vector<double> uKnots(boundaries[0]->pole_num() - boundaries[0]->degree() - 1);
    std::copy(boundaries[0]->knots().begin() + boundaries[0]->degree() + 1,
              boundaries[0]->knots().begin() + boundaries[0]->pole_num(), uKnots.begin());
    T->refine_u_knot(std::move(uKnots));

    std::vector<double> vKnots(boundaries[1]->pole_num() - boundaries[1]->degree() - 1);
    std::copy(boundaries[1]->knots().begin() + boundaries[1]->degree() + 1,
              boundaries[1]->knots().begin() + boundaries[1]->pole_num(), vKnots.begin());
    T->refine_v_knot(std::move(vKnots));

    // S = S1 + S2 - T
    for (std::size_t i = 0; i < boundaries[0]->pole_num(); i++)
        for (std::size_t j = 0; j < boundaries[1]->pole_num(); j++)
            dstPoles[i * boundaries[1]->pole_num() + j] -= T->pole(i, j);

    auto S = std::make_shared<Geom_BSplineSurface<3>>(std::move(dstPoles), T->move_u_knots(), T->move_v_knots(),
                                                      T->u_degree(), T->v_degree());
    return S;
}

std::shared_ptr<Geom_BSplineSurface<3>> Alg_FairingEnergy_Covering::cover(
    std::vector<std::shared_ptr<Geom_BSplineCurve<3>>> &&curves)
{
    auto boundaries = compat_curves(fix_curves_order(std::move(curves)), true);
    if (boundaries.size() == 0)
        return nullptr;

    if (boundaries.size() >= 4)
    {
        Alg_BSplineSurfaceEnergy_Approximator<3> approximator;
        approximator.uDegree = uDegree;
        approximator.vDegree = vDegree;
        approximator.wStretch = wStretch;
        approximator.wBending = wBending;

        std::size_t samples = boundarySampleNum;
        for (std::size_t i = 0; i < boundaries.size(); i++)
        {
            auto curve = boundaries[i];
            auto first = curve->first_param();
            auto last = curve->last_param();
            auto step = (last - first) / (samples - 1);
            for (std::size_t j = 0; j < samples; j++)
            {
                auto param = first + j * step;
                auto point = curve->value(param);
                approximator.data.push_back(point);
            }
        }

        auto result = approximator.approximate(uPoleNum, vPoleNum);
        params = std::move(approximator.params);
        data = std::move(approximator.data);
        return result;
    }

    auto p = boundaries[0]->degree();
    auto q = boundaries[1]->degree();

    // 4 个基函数相乘，要确保精确，至少需要 2 * p 阶
    static Integrator_Gauss integrator(p + q);

    auto uKnots = boundaries[0]->knots();
    auto vKnots = boundaries[1]->knots();
    auto n1 = boundaries[0]->knot_num() - p - 2;
    auto n2 = boundaries[1]->knot_num() - q - 2;

    // 计算基函数的向量
    Eigen::MatrixXd Bu0(n1 + 1, n1 + 1), Bu1(n1 + 1, n1 + 1), Bu2(n1 + 1, n1 + 1);
    Eigen::MatrixXd Bv0(n2 + 1, n2 + 1), Bv1(n2 + 1, n2 + 1), Bv2(n2 + 1, n2 + 1);
    Bu0.setZero();
    Bu1.setZero();
    Bu2.setZero();
    Bv0.setZero();
    Bv1.setZero();
    Bv2.setZero();

    auto &weights = integrator.weights;
    auto &nodes = integrator.nodes;
    for (std::size_t i = p; i < n1 + 1; i++)
    {
        double mid = (uKnots[i + 1] + uKnots[i]) / 2;
        double ratio = (uKnots[i + 1] - uKnots[i]) / 2;

        // 防御重节点，此时区间为零，不用积分
        if (std::abs(ratio) > std::numeric_limits<double>::epsilon())
        {
            for (std::size_t j = 0; j < nodes.size(); j++)
            {
                auto derivs = deriv_basis_funcs(uKnots.data(), p, i, 2, mid + nodes[j] * ratio);
                for (std::size_t ii = 0; ii < p + 1; ii++)
                {
                    for (std::size_t jj = 0; jj < p + 1; jj++)
                    {
                        Bu0(i - p + ii, i - p + jj) += derivs[0][ii] * derivs[0][jj] * weights[j];
                        Bu1(i - p + ii, i - p + jj) += derivs[1][ii] * derivs[1][jj] * weights[j];
                        Bu2(i - p + ii, i - p + jj) += derivs[2][ii] * derivs[2][jj] * weights[j];
                    }
                }
            }
        }
    }

    for (std::size_t i = q; i < n2 + 1; i++)
    {
        double mid = (vKnots[i + 1] + vKnots[i]) / 2;
        double ratio = (vKnots[i + 1] - vKnots[i]) / 2;

        // 防御重节点，此时区间为零，不用积分
        if (std::abs(ratio) > std::numeric_limits<double>::epsilon())
        {
            for (std::size_t j = 0; j < nodes.size(); j++)
            {
                auto derivs = deriv_basis_funcs(vKnots.data(), q, i, 2, mid + nodes[j] * ratio);
                for (std::size_t ii = 0; ii < q + 1; ii++)
                {
                    for (std::size_t jj = 0; jj < q + 1; jj++)
                    {
                        Bv0(i - q + ii, i - q + jj) += derivs[0][ii] * derivs[0][jj] * weights[j];
                        Bv1(i - q + ii, i - q + jj) += derivs[1][ii] * derivs[1][jj] * weights[j];
                        Bv2(i - q + ii, i - q + jj) += derivs[2][ii] * derivs[2][jj] * weights[j];
                    }
                }
            }
        }
    }

    auto nu = boundaries[0]->pole_num();
    auto nv = boundaries[1]->pole_num();

    auto N = nu * nv;
    auto Nint = (nu - 2) * (nv - 2);

    // 构建对应关系
    auto m2v = [&](int i, int j) -> int {
        // 边界情况
        if (i == 0)
            return Nint + j;
        if (i == nu - 1)
            return Nint + nv + j;
        if (j == 0)
            return Nint + nv * 2 + i - 1;
        if (j == nv - 1)
            return Nint + nv * 2 + nu - 2 + i - 1;

        return (i - 1) * (nv - 2) + j - 1;
    };

    auto v2m = [&](int ind, int &i, int &j) {
        // 内部情况
        if (ind < Nint)
        {
            i = ind / (nv - 2) + 1;
            j = ind % (nv - 2) + 1;
        }
        // 边界情况
        else if (ind < Nint + nv)
        {
            i = 0;
            j = ind - Nint;
        }
        else if (ind < Nint + nv * 2)
        {
            i = nu - 1;
            j = ind - (Nint + nv);
        }
        else if (ind < Nint + nv * 2 + nu - 2)
        {
            i = ind - (Nint + nv * 2) + 1;
            j = 0;
        }
        else
        {
            i = ind - (Nint + nv * 2 + nu - 2) + 1;
            j = nv - 1;
        }
    };

    std::vector<Eigen::Triplet<double>> tripletsA;
    std::vector<Eigen::Triplet<double>> tripletsB;

    // AP1 = BQ
    for (int i = 0; i < N; i++)
    {
        int ui, vi;
        v2m(i, ui, vi);

        for (int j = 0; j < N; j++)
        {
            int uj, vj;
            v2m(j, uj, vj);

            // 能量项
            double Euu = Bu2(ui, uj) * Bv0(vi, vj);
            double Euv = Bu1(ui, uj) * Bv1(vi, vj);
            double Evv = Bu0(ui, uj) * Bv2(vi, vj);

            double r = Euu + 2 * Euv + Evv;
            if (j < Nint)
                tripletsA.emplace_back(i, j, r);
            else
                tripletsB.emplace_back(i, j - Nint, -r);
        }
    }

    Eigen::SparseMatrix<double> A(N, Nint);
    A.setFromTriplets(tripletsA.begin(), tripletsA.end());

    Eigen::SparseMatrix<double> B(N, N - Nint);
    B.setFromTriplets(tripletsB.begin(), tripletsB.end());

    // 填充边界控制点
    Eigen::MatrixXd Q(N - Nint, 3);
    std::vector<Vec3> dstPoles(nu * nv);
    for (std::size_t i = 0; i < nu; i++)
    {
        auto P1 = boundaries[0]->pole(i);
        auto P2 = boundaries[2]->pole(boundaries[2]->pole_num() - 1 - i);

        auto ind1 = m2v(i, 0) - Nint;
        auto ind2 = m2v(i, nv - 1) - Nint;

        Q(ind1, 0) = P1.x();
        Q(ind1, 1) = P1.y();
        Q(ind1, 2) = P1.z();

        Q(ind2, 0) = P2.x();
        Q(ind2, 1) = P2.y();
        Q(ind2, 2) = P2.z();

        dstPoles[i * nv] = P1;
        dstPoles[(i + 1) * nv - 1] = P2;
    }

    for (std::size_t i = 0; i < nv; i++)
    {
        auto P1 = boundaries[3]->pole(boundaries[3]->pole_num() - 1 - i);
        auto P2 = boundaries[1]->pole(i);

        auto ind1 = m2v(0, i) - Nint;
        auto ind2 = m2v(nu - 1, i) - Nint;

        Q(ind1, 0) = P1.x();
        Q(ind1, 1) = P1.y();
        Q(ind1, 2) = P1.z();

        Q(ind2, 0) = P2.x();
        Q(ind2, 1) = P2.y();
        Q(ind2, 2) = P2.z();

        dstPoles[i] = P1;
        dstPoles[(nu - 1) * nv + i] = P2;
    }

    // 稀疏 QR 分解
    Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
    solver.compute(A);

    Eigen::MatrixXd x = solver.solve(B * Q);
    for (int k = 0; k < Nint; k++)
    {
        int i, j;
        v2m(k, i, j);

        dstPoles[i * nv + j] = Vec3(x(k, 0), x(k, 1), x(k, 2));
    }

    auto S = std::make_shared<Geom_BSplineSurface<3>>(std::move(dstPoles), std::move(uKnots), std::move(vKnots), p, q);
    return S;
}

std::shared_ptr<Geom_BSplineSurface<3>> Alg_FairingPIA_Covering::cover(
    std::vector<std::shared_ptr<Geom_BSplineCurve<3>>> &&curves)
{
    std::shared_ptr<Geom_BSplineSurface<3>> result;
    if (curves.size() <= 4)
    {
        Alg_Coons_Covering coons;
        auto baseSurface = coons.cover(std::move(curves));
        if (baseSurface == nullptr)
            return nullptr;

        auto [U1, U2, V1, V2] = baseSurface->bounds();

        std::size_t samples = boundarySampleNum;
        Alg_Surface_FairingPIA<3> fairingPIA;
        fairingPIA.surface = baseSurface;
        fairingPIA.data.resize(samples * 4);
        fairingPIA.params.resize(samples * 4);

        auto stepU = (U2 - U1) / (samples - 1);
        auto stepV = (V2 - V1) / (samples - 1);
        for (std::size_t i = 0; i < samples; i++)
        {
            fairingPIA.data[i] = baseSurface->value(U1 + i * stepU, V1);
            fairingPIA.data[i + samples] = baseSurface->value(U1 + i * stepU, V2);
            fairingPIA.data[i + 2 * samples] = baseSurface->value(U1, V1 + i * stepV);
            fairingPIA.data[i + 3 * samples] = baseSurface->value(U2, V1 + i * stepV);
            fairingPIA.params[i] = Vec2{U1 + i * stepU, V1};
            fairingPIA.params[i + samples] = Vec2{U1 + i * stepU, V2};
            fairingPIA.params[i + 2 * samples] = Vec2{U1, V1 + i * stepV};
            fairingPIA.params[i + 3 * samples] = Vec2{U2, V1 + i * stepV};
        }

        auto uRange = baseSurface->u_degree() + 1;
        auto vRange = baseSurface->v_degree() + 1;
        auto un = baseSurface->u_pole_num();
        auto vn = baseSurface->v_pole_num();

        fairingPIA.wBending.resize(un * vn, wBendingBoundary);
        fairingPIA.wStretch.resize(un * vn, wStretchBoundary);
        for (std::size_t i = 0; i < un; i++)
        {
            for (std::size_t j = 0; j < vn; j++)
            {
                if (i < uRange || i > un - uRange - 1 || j < vRange || j > vn - vRange - 1)
                    continue;

                fairingPIA.wBending[i * vn + j] = wBendingInterior;
                fairingPIA.wStretch[i * vn + j] = wStretchInterior;
            }
        }

        auto iter = fairingPIA.generate();
        for (auto &sur : iter(iterNum))
            result = sur;
    }
    else
    {
        Alg_DBS_Approximator<3> approximator;
        approximator.uPoleNum = uPoleNum;
        approximator.vPoleNum = vPoleNum;
        approximator.uDegree = uDegree;
        approximator.vDegree = vDegree;
        approximator.wStretch = 1e-4;
        approximator.wBending = 1e-3;

        std::size_t samples = boundarySampleNum;
        for (std::size_t i = 0; i < curves.size(); i++)
        {
            auto curve = curves[i];
            auto first = curve->first_param();
            auto last = curve->last_param();
            auto step = (last - first) / (samples - 1);
            for (std::size_t j = 0; j < samples; j++)
            {
                auto param = first + j * step;
                auto point = curve->value(param);
                approximator.data.push_back(point);
            }
        }

        // 使用动态基曲面构造初始曲面
        std::shared_ptr<Geom_BSplineSurface<3>> baseSurface;
        auto iterBase = approximator.generate();
        for (auto &sur : iterBase(3))
            baseSurface = sur;

        auto un = baseSurface->u_pole_num();
        auto vn = baseSurface->v_pole_num();
        auto uDegree = baseSurface->u_degree();
        auto vDegree = baseSurface->v_degree();
        auto &uKnots = baseSurface->u_knots();
        auto &vKnots = baseSurface->v_knots();

        Alg_Surface_FairingPIA<3> fairingPIA;
        fairingPIA.surface = baseSurface;
        fairingPIA.data = std::move(approximator.data);
        fairingPIA.params = std::move(approximator.params);
        fairingPIA.wBending.resize(un * vn, wBendingInterior);
        fairingPIA.wStretch.resize(un * vn, wStretchInterior);

        // 投影采样点，考虑受影响的控制点权
        for (std::size_t i = 0; i < fairingPIA.params.size(); i++)
        {
            auto P = fairingPIA.params[i];
            auto uSpan = find_span(uKnots.data(), uKnots.size(), uDegree, P.x());
            auto vSpan = find_span(vKnots.data(), vKnots.size(), vDegree, P.y());

            for (std::size_t j = uSpan - uDegree; j < uSpan + 1; j++)
            {
                for (std::size_t k = vSpan - vDegree; k < vSpan + 1; k++)
                {
                    fairingPIA.wBending[j * vn + k] = wBendingBoundary;
                    fairingPIA.wStretch[j * vn + k] = wStretchBoundary;
                }
            }
        }

        auto iterResult = fairingPIA.generate();
        for (auto &sur : iterResult(iterNum))
            result = sur;

        params = std::move(fairingPIA.params);
        data = std::move(fairingPIA.data);
    }

    return result;
}